function [TMatrix,EPara] = MHMMESTIMATE(Dyn,Values,varargin)

pseudoEParacounts = false;
pseudoTMatrixcounts = false;
TMatrix = [];
EPara = [];
DynLen = length(Dyn);
if length(Values) ~= DynLen
    error('stats:hmmestimate:InputSizeMismatch',...
          'Input sequence and states must be the same length.');
end

uniqueValues = unique(Values);
numValues = length(uniqueValues);
States = zeros(DynLen,1);
for count = 1:numValues
    States(Values==(uniqueValues(count)))=count;
end

numStates = numValues;
numEPara = 2; % mean % sigma

TMatrix = zeros(numStates);
EPara = zeros(numStates, numEPara);

% count up the transitions from the state path
for count = 1:DynLen-1
    TMatrix(States(count),States(count+1)) = TMatrix(States(count),States(count+1)) + 1;
end

% and calculate the emissions parameters for each state
for count = 1:numStates
    Group = Dyn(States == count);
    EPara(count, 1) = mean(Group); %mean
    EPara(count, 2) = std(Group); %sigma
end

trRowSum = sum(TMatrix,2);
% if we don't have any values then report zeros instead of NaNs.
trRowSum(trRowSum == 0) = -inf;

% normalize to give frequency estimate.
TMatrix = TMatrix./repmat(trRowSum,1,numStates);
